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ABSTRACT 



Two methods of obtaining sensitivity data were simulated on an 
electronic computer for the purpose of comparing the accuracy of the 
estimates of the parameters of an underlying cumulative normal response 
function. The first method simulated the standard Bruceton procedure 
while the second used a modified binary search routine with a portion of 
the sample in order to obtain maximum likelihood estimates of the input 
parameters for use in a follow-on Bruceton test. 

The results showed both methods to be effective in estimating the 
mean but with slightly more variability in the estimates obtained by the 
second procedure. Both methods underestimated the standard deviation - 
again with more variability in the estimates obtained by the second 
procedure. When the prior parameter estimates were unknown and the 
applicable stimulus level bounded, the second method yielded estimates 
favorably comparable to those expected from the Bruceton procedure with 
suitable prior input estimates. 
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I. INTRODUCTION 



Frequently a statistician is faced with the problem of determining 
the level of a stimulus which critically affects the performance of a 
device. The nature of the testing to be discussed is such that once 
some positive level of the stimulus is applied to the device either a 
response or a non-response can be immediately observed and, in either 
case, the device is altered so that a bonafide result cannot be obtained 
from a second test. Tests of this type are known as sensitivity tests. 

One of the many problems besetting those involved in explosives 
research is that of providing measures and specifying rules to provide 
for the safe handling and transportation of explosives. Many different 
types of sensitivity testing apparatus have been developed for laboratory 
use, the most common being those that subject some quantity of explosive 
to the impact load of a falling drop-weight from some controllable height. 
At least as late as October 1965 there remained two important physical 
problems to be solved; namely, that of establishing a measure of stimulus 
not highly apparatus -dependent and then that of translation of these 
results to safe handling rules [1]. These problems are not addressed in 
this paper but should be kept in mind when considering the overall problem. 

In the early 1940’s, a technique for obtaining sensitivity data was 
developed and used in explosives research at the Explosives Research 
Laboratory, Bruceton, Pennsylvania which has come to be called synony- 
mously, the Bruceton, Staircase, or ^*Up and Down” Method. 

The aim of this method of testing is to increase the accuracy with 
which certain critical values of the stimulus may be estimated, notably 
the median (or mean) and standard deviation. The accuracy of the method 
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depends in part on the stimulus level at which the first item is tested 
and the interval spacing for subsequent levels of testing [2]. 

VJhen the stimulus levels mentioned above cannot be determined prior 
to testing or when little confidence is placed on the available estimates, 
a preliminary (or search) phase of testing may be desirable to obtain 
maximum likelihood estimates prior to employing the Bruceton Method with 
the remainder of the sample, A procedure to do this is offered as an 
alternative method . 

The comparative accuracies of the two techniques were examined 
through the use of simulation conducted on a high-speed electronic 
computer. All parameters and estimates considered as inputs to the 
simulation were kept within ranges for which the Bruceton Method is 
considered to yield accurate results [2], 
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II. THE MODEL 



Let X be an applied stimulus level (x 0 [o,co)) and y = y(x) be the 
associated response (ye^o>l^ where **o” denotes no response and ” 1 *' 
denotes response). At any given stimulus level consider y to be the 
realization of a Bernoulli random variable, Y, with response 
probability 

p(x) = Prob (Y = l|x) 

The function p(x) is called the response function and is further 
specified as 

p(x) = 0 xe[o,a] 

0 < p(x) < 1 xe(a,b) 

and p(x) 1 xe[b,oo) 

The intervals [o,a], (a,b) and [b,co) are called the zero-response 
region, the mixed-response region, and the one-response region 
respectively. It is assumed that p(x) is a monotonely increasing 
function for stimulus values in the mixed-response region. Thus, 
p(x) can be considered as the cumulative distribution function for a 
random variable X such that 

p(x) = Prob (X < x). [3] 

In this context the random variable X can be interpreted as a thres- 
hold stimulus level, thus 

Prob (Y = ijx) = Prob (X < x) = p(x) 

and Prob (Y = o]x) = Prob (X > x) = 1 - p(x), [3] 

2 

It is assumed the X is distributed Normal (|i,o ); that is 

p(x) = cp(x|ii,a ) 

2 

where cp(x|(i,a ) represents the cumulative normal distribution with mean 
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(a and variance o . 



In particular 



Prob (x < ij) = p(^) = 0.5. [3] 
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III. TESTING METHODS 



A, BRUCETON METHOD 
1 . Description 

Based on intuition or past experiments, the experimenter selects 
a priori estimates of \l and a. Call these estimates and and let 



d = 



The 



experimenter tests the first item at or near |j.^ . If there 



is a response the second item is tested at a level d units below 
otherwise the second item is tested at a level d units above In 

the same manner, each of the remaining items is tested at a level d 
units above or below the previous test level according as there was not 
or there was a response observed for the previous test. Thus the 
sample is concentrated about the mean and one would expect nearly equal 
numbers of responses and non-responses. In fact, the number of non- 
responses at any level will not differ by more than one from the 
number of responses at the next higher level [2]. 

Let N denote the total number of observations of the less 
frequent event and n^ ,n^^ ,n 2 , • • *nj^ denote the frequencies of this event 
at each level where n^ corresponds to the lowest level and n^ the 
highest level at which the less frequent event occurs. 

The final estimates of fi and a are based on the first two 
moments of the stimulus levels. Since the intervals are equally 
spaced, these moments can be computed in terms of the sums 

A = S i n. 

. 1 
1 



and 



B = S i^ n . . 
i 
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Let \ji be the estimate of (j. by this method. Then 

2 = x' + d(|±£) 

where x* represents the lowest level at which the less frequent event 
occurs [2]. The plus sign is used when the analysis is based on non- 
responses, and the minus sign when it is based on responses [2]. 

2 2 

If (NB-A )/N > .3 the sample standard deviation is 

s = 1.620 d -■ + . 029 ') 

Otherwise, a more elaborate calculation must be employed and is 
described in Ref. 2. 

To obtain confidence intervals, estimates of the standard 
deviations of the sample mean and sample standard deviation, say s^ 
and respectively, are given by 

Gs 

and 

Hs 

"" ~7n 

where the factors G and H are dependent on the ratio “ and the 
position of the mean relative to the testing levels. Plots of these 
factors are available in Ref. 1. 

2 . Discussion 

Only rarely is the threshold stimulus Z normally distributed. 
It is usually the case that some scale transformation of Z, say X, is 
made so that X is normally distributed in the vicinity of the mean. 
This transformation is done prior to testing to determine and a^. 
Only after all analysis is completed are the values scaled back to the 
original stimulus measure [2]. 
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The size of the sample is critical to the accuracy of the 
estimation. Note that at most only half of the sample is used in the 
analysis so that, for example, if thirty items are tested the maximum 
possible value of N is fifteen. The analysis is based on large sample 
theory which in the case mentioned would be applied to a sample of size 
fifteen [2 ] [4] . 

Unless normality of the variate is assured this method does 
not yield accurate results for the small and large percentage points. 
This is unfortunate since in most applications one would be more 
interested in a small percentage point as a measure of safety and a 
large percentage point as a measure of reliability. At any rate, an 
estimate of a percentage point j is 

A 

j = + ks 

where k is chosen from tables of the standard normal deviate to give 
the desired percentage [2]. One could then conduct tests in the 
vicinity of this value to refine the estimate. 

B. BRUCETON METHOD PRECEDED BY SEARCH 
1 . Description 

In the event that a priori estimates of ^ and a are not 
available some economic method of attaining these estimates is desired. 
A method proposed and described below is a modified binary search 
technique. 

Again, the assumption is that the threshold stimulus (or 
some transformation of it) is normally distributed and p(x) can be 
represented by a cumulative normal distribution. 

As noted from the model 

Prob (Y = 0|x < a) = 1 
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and 



Prob (Y = l|x>b) = 1. 

The first step in the procedure, then, is to select values for a and b, 
(In the case of complete uncertainty these could be the limiting values 
of the testing apparatus) and commence the binary search starting at 

x= (a + b)/2. 

If p(x) were a step function, repetition of this method would locate 
the step in an interval of any desired length. In general, however, 
the mixed-response region has non-zero width and a non-response would 
merely indicate that the applied stimulus is in the mixed response 

region or below while a response would indicate that it was in the 

mixed response region or above. 

If a test at x^^ yields a response and a test at x^ yields a 

non-response while < x^ it is certain that both x^^ and x^ are in the 

mixed response region. This condition is called a response inversion 
and is the basic indicator for the modified binary search technique. 

The description of the procedure is best followed by referring to 
Figures 1 through 4. 

Sequence S* is a cyclic one indicating that a reduction in step 
size should be taken. Test levels are selected attempting to reproduce 
this sequence. Failure to do this results in the basic inversion 
sequence S^. Tests are then made at the end of this sequence to result 
in one of three ternninal situations S^, S 2 > or S^. In the event the 
mixed response region is relatively narrow and near a or b , several 
binary reductions may be necessary to reproduce S* or one of the 
terminal situations. These circumstances are represented by and 
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STARTING SEQUENCE 
Figure 1 
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S^J SEQUENCE 
Figure 2 
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S* SEQUENCE 
Figure 3 
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S SEQUENCE 

Lf 

Figure 4 
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Maximum likelihood estimates of \i and a are available for 



sequences and and developed as described below [ 3 ]. 

2 . Discussion 

It is assumed that all trials are independent. Thus the 
probability of the sequence is 

Prob (Sp = Prob (Y^=0, ¥^=0, ¥ 3 =!, Y^=0, ¥ 3 =!, Y^=llx^,X 2 ,X 3 ,X^,X 3 ,X^) 

6 

= rr Prob (Y. = y . lx . ) 

. T 1 1 ' 1 

1 = 1 



where 

Prob (Y. = y.lx.) = cp(x.) if y. =1 

= 1 - cp(x^) if = 0 

and « 

cp(x.) = Prob (X. < X.) = f — — e ^ dx. 

" J 



Maximum likelihood estimates for i_i and o can then be established 
using standard normal tables for each of the terminal situations. 

These estimates are indicated on Figure 5. 
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RESPONSE SEQUENCES AND PARAMETER ESTIMATES 
Figure 5 



20 



IV. SIMUIAI.ON 



A. DESCRIPTION 

All simulated experiments were conducted on an IBM 360/67 computer 
using the FORTRAN IV programming language. The basic program is 
attached. The response function p(x) used was cumulative normal v;ith 
ji = 30 and a = 3 . 

The sample size was kept at seventy for each experiment to provide 
some assurance that the analytical sample would be suitable for large 
sample analysis. 

The basic test procedure was to draw a random number on the unit 
interval and compare this to F(x), a function of a standard normal 
variate specified as 



and 



where 




if X < 0 



if X > 0 



(The function subprogram erf is an IBM-supplied subprogram.) If the 

random number was less than or equal to F(x^) then a response was 
t h 

counted for the i level; otherwise a non-response was counted. 

Six different cases were tested using the straight Bruceton pro- 
cedure (METHOD 1) with two different input estimates of and three 
different input estimates of a. Case 1 considered exact estimates; 
i.e., = |J and = a. Case 2 considered = ^i-6 and = a* 
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Cases 3 and 4 considered and = a/2, 2a respectively while 

Cases 5 and 6 repeated Cases 3 and 4 except |a^ = (J-6. For each of the 
six cases 1000 experiments were conducted each utilizing a different 
sequence of random numbers . 

The search procedure (METHOD 2) was then incorporated into each 
of the above six cases using the a prior estimates, \j.^ and a^j to 
determine estimates for stimulus levels a and b and thereby the size 
of the binary reduction as indicated in Figure 1, The program then 
followed the flow shown in Figures 1 through 4 until either a terminal 
sequence was reached or the search was arbitrarily terminated as dis- 
cussed in subparagraph C below. The Bruceton procedure was then used 
until the sample was exhausted. 

The final case. Case 7, indicated complete lack of knowledge of ^ 
and a but considered the upper and lower stimulus level limits of the 
test apparatus to be 100 and 0 respectively. 

B. MEASURES OF EFFECTIVENESS 

At the completion of all experiments for each case, several 
measures were obtained for comparison. First, average values of the 
parameters were determined to be 

ji = E C./N 

i 

and 

a = E a./N 
i 

A A th 

where and a^ are the a posteriori estimates of |J and a for the i 

experiment and n, the number of experiments used. Next, as measures of 

variability 
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and 



s^A = E (m. - ^)^/n-l 
^ i 

2 ^ 

s A = E (o^ - o)/n-l 

were calculated. In addition, the program listed the maximum and 
minimum estimates of both |j and a. 

C. DISCUSSION 

In Chapter III it was noted that sequences S*, S , and S are 

U Li 

cyclic. In order to simplify the program it was necessary to 
artificially terminate these situations at some point and calculate 
the input values for the Bruceton test. The estimate of (j. used was 
Mg 

where and X 2 are adjacent testing levels and X 2 > Xj^ with y^^ = 0 
and y 2 == 1. The estimate of a used was 

for Cases 1 through 6 and 




for Case 7. The former estimate of a was chosen arbitrarily while the 
latter estimate was based on the estimate of the mixed response region 
being 6a. While the number of terminations of this type was insignifi 
cant for the first six search cases, in the final case over 600 experi 
ments were thus terminated requiring the program to be expanded to 
permit more recycling. The point is that the artificial termination 
does not represent the search procedure. This problem would not arise 
in field experimentation until either the sample was exhausted or the 
step size reduction of stimulus level indicated was too narrow to be 
measured or controlled by the test apparatus. 
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Also in the interest of program simplification those experiments 



for which 




< .3 



were not used for analyses. This limitation invalidated the measures 
of effectiveness for the Bruceton cases where = 2 q. 



D. RESULTS 

The results of the simulation are listed in Table I. It is 
questionable that the measures listed under Method 1 are valid for 
Cases 4 and 6 in that only ,381 and .393 of the possible experiments 
were used. These two cases and Case 4 under Method 2 (where .661 of 
the possible experiments were used) are the only ones for which 
a > a. 

In general the extreme estimates are more widely separated and the 
variability of a is greater in Method 2. 

Estimates of |i range from 27.8823 to 31.7647 for Method 1 and 
27.937 to 31.91 for Method 2. 

Estimates of a range from .8741 to 6.5027 for Method 1 and .3498 
to 9.8328 for Method 2. 

The lowest average (i, 29.9113, was obtained under Method 1, Case 5, 
while the highest average ja, 30.1175, was obtained under Method 2, 

Case 3 . 

The lowest average a, 2.3748, was obtained under Method 2, Case 5, 
while the highest average a, 2.9474, was obtained under Method 1, 

Case 5. (Case 6 is not counted under Method 1 nor is Case 4 under 
both methods . ) 
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TABLE OF EXPERIMENTAL RESULTS 





METHOD 1 


METHOD 2 


A 

H 


A 

a 


A 

H 


A 

a 


CASE 1 










Hi = 30 


AVE 


30.0067 


2.8320 


30.0117 


2.8609 


Oi =3 


MAX 


31.7647 


5.7904 


31.7813 


5.9343 


a = 18 


MIN 


28.5000 


1.6089 


28.2187 


1. 1241 


b = 42 


VAR 


.2523 


.4128 


.2514 


.5831 


CASE 2 










Hi = 24 


AVE 


29.9641 


2.9040 


30.0317 


2.8819 


Oi = 3 


MAX 


31.6765 


5.8249 


31.6875 


9.1369 


a = 12 


MIN 


28.3235 


1.6250 


28.1976 


.9512 


b = 36 


VAR 


.2656 


.4225 


.2666 


.6336 


CASE 3 










Hi = 30 


AVE 


30.0295 


2.7216 


30.1175 


2.8615 


Oi = 1.5 


MAX 


31.6071 


6.1197 


31.9100 


7.7997 


a = 24 


MIN 


28.4118 


.8741 


28.5950 


.8697 


b = 36 


; VAR 


.2046 


.7409" 


.2693 


.9081 


CASE 4 










Hi = 30 


AVE 


29.9683 


3 . 5424 


29.9750 


3.0721 


Qj = 6 


MAX 


31.4571 


6.0266 


31.6875 


6.3569 


a = 6 


MIN 


28.0286 


-- 


27.9370 


1.6170 


b = 54 


VAR 


.2574 


.4639 


.2619 


.4522 


CASE 5 










Hi = 24 


AVE 


29.9113 


2.9474 


29.9363 


2.3748 


Oi = 1.5 


MAX 


31.4773 


5.9257 


31.4063 


7.9507 


a = 18 


MIN 


28.2353 


.9452 


28.1961 


.3498 


b = 30 


VAR 


.2220 


.8748 


.2184 


1.4889 


CASE 6 










Hi = 24 


AVE 


29.9493 


3.5438 


30.0247 


2.8398 


a = 6 


MAX 


31.4118 


6.5027 


31.5756 


6.4201 


a = 0 


MIN 


27.8823 


-- 


28.2552 


1.1252 


00 

II 


VAR 


.2639 


.4785 


.2490 


.6300 


CASE 7 










-- 


AVE 


-- 


-- 


30.0123 


2.7280 


-- 


MAX 


-- 


-- 


31.8229 


9.8328 


a = 0 


MIN 


-- 


-- 


27.9541 


.5082 


b = 100 


VAR 


-- 


-- 


.2628 


2.1041 



TABLE I 
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V. CONCLUSIONS AND RECOMMENDATIONS 



A. CONCLUSIONS 

1. Estimation of the Mean 

Both methods estimate the mean effectively. 

2 . Estimation of the Standard Deviation 

Both methods tend to under-estimate the standard deviation with 
no predictable bias and are therefore unsuitable for use in safety or 
reliability statements. This conclusion agrees with the findings of 
Hampton [4] as it pertains to the Bruceton Method. 

3 . Extension of the Search Phase for the Starting Sequence 

Termination of the search phase with sequence S^^ in the starting 

sequence (see Figure 1) may yield estimates of a greater than twice 
the actual value. To avoid this it is advisable to extend the search 
phase as described in Ref. 3, 

4. Use of Search Technique 

The search procedure should be used in those cases where there 
is not independent evidence that the estimate of a is within the range 
for which the Bruceton Method is recommended (i.e., a/2 < < 2 q) . 

B. RECOMMENDATIONS 

Further testing of Method 2 is recommended under the circumstances 
listed below. 

1 • Reduction of Sample Size 

It would be of in‘*erest to reduce the sample size to the point 
where the effective sample is small, say 15, and compare the Bruceton 
procedure with the search procedure using the entire sample for the 
search . 
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2 . Random Selection of Response Function Parameters 



A more valid test of both methods would be achieved by 
randomally selecting values of |j and a over some range and using 
the on-line computer facility to conduct the simulation. 
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COMPUTER PR^GR^M 



THIS program SI»^"ULATES SENSITIVITY TFSTING BY BOTfl THE BRllC- 
FTON method (WHEN IANY=0) AND THE BRUCETON METHOD PPPCEDED 
BY THE MODIFTED BINARY SEARCH (WHEN IANY=1).THE UNDEPl.YIMG 
RESPONSE FUNTION IS CUMULATIVE NORMAL {BO,B).THr- IHDijT EST- 
IMATES OF THE MEAN AND THE STANDARD DEVIATION ARE CALLED 
EXMU AND EXSIG RESPECTIVELY. 

THE PRINCIPLE VAPIARLF NAMES ARE AS FOLLOWS... 

AACT IS THC STIMULUS VALUE AT THE UPPER LIMIT 0^ THE 
MIXED RESPONSE »EGION. 

BACT IS THF STIMULUS VALUE AT THE LOWER LIMIT OF THE 
MIXED RESPONSE REGION. 

A AND B ARE ESTIMATFS OF AACT AND RACT RES PECT I V FLY . 

X(J) IS THE STIMULUS LEVEL OF THE JTH. STIMULUS. 

IXO(J) IS the cumulative count of NON-RESPONSFS at X(J). 

IXX(J) IS the CUMUIATIVF COUNT OF RESPONSES AT X(J) 

IS IS THE sample SIZE. 

NU IS THE ENTRY NUMBER FOR THE RANDOM NUMBER GENFRAT0R» 

UN I F . 

N COUNTS THE NUMBER OF FXPERIMENTS. 

RN IS THE RANDOM NUMBER ON (0,1) RETURNFO BY UNIF. 

FOFX IS THE VALUE OF THE RESPONSE FUNCTION RETURNED BY 
SUBPROGRAMS XNCOF AND SNCDF. 

ISUMO IS THE TOTAL NUMBER OF NON-RESPONSES FOR ONE EXPER- 
IMENT. 

ISHMX IS THE TOTAL NUMBER OF RESPONSES FOR ONE EXPERIMENT 

NT IS THE MINIMUM OF ISUMO AND ISUMX. 

NS(J) IS the FPEOUENCY DF THE LESS FPEOUENT EVENT AT X(J) 

NG(J) REARRANGES NS(J) SO THAT NGm=MS(M WHEPT; X(I) IS 
THE LOWEST STIMULUS LEVEL AT WHICH THE LESS FREQUENT 
EVENT OCCURS. 

AR(J) IS USED TO CALCULATE THE first MOMENT ,SUMAP . 

BR(J) IS USED TO CALCULATE THE SECOND MOMENT , SUM'^R . 

YPRIME IS THE LOWEST LEVEL AT WHICH THE LESS FREQUENT 
EVENT OCCURS 

XMUEST IS THE FINAL ESTIMATF 0= THE TRUE M-AN,XMU. 

DEVEST IS THE FINAL ESTIMATE OF THE TRUE STANDARD DEVIAT- 
lONtXSIG. 

EMU(J) IS THF difference OF XMUEST AND XMU. 

EDEV(J) IS THE D»FFERENCE OF DEVEST AND XSIG. 

SAMAVM AND SAMAVD ARE THE SAMPLE AVERAGE ERRORS OF XMUEST 
AND DEVEST RESPECTI VFLY. 

SAMSQM AND SAMSQO ARE THE AVERAGE Mfan SQUARE frrprs OF 
XMUEST AND DEVEST RESPECTIVELY. 

NOGO IS THE NUMBER OF EXPERIMENTS NOT USED IN THE FINAL 
ANALYS IS 



dimension arrays and format 

SIMULATE BRUCETON FIRST THEN SEARCH 

IANY=0 

69 THING=0. 

INITIALIZE INTERNAL AND OUTPUT VARIABLES 

SET EXPERIMENT COUNTER, SAMPLE SIZE COUNTfr, aMO NU. 

LC0UNT=1 
NU=l??71 
103 N=1 

IF( lANY.EQ.O) NBR = 0 
IF( lANY.EO.l ) NBR = 1 

SET INPUT VARIABLES 

XMU=30. 

XSTG=3 . 

A=0. 

B=100. 

EXMU=50. 
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EXSir,= 12. 5 
A= FXMU-IO’f'PXSIG 
B= FXMU + IO*f=XSIG 
Xl=(A+B) /2. 

INC = 0 



PROVIDE BRANCH TO STANDARD BRUCETON 
IF(NBR.EO.O) GO TO 33 

CONDUCT SEARCH 



CALL UNIF(RN,NU) 
FOFX=XNCDF(Xl tXMUtXSIG) 

IF( RN.LE.FO«^X ) GO TO 9500 
X2=( B+Xl) /2. 

NBR=NBR+1 
CALL UNIP(RN,N'J) 
FOFX=XNCDF( X?,XMU, XSIG> 
IF{RN.LE.FOPX) GO TO 9250 
X3=( R+X2) /2. 

NBR=NBR+1 

CALL UNIF(RN,NU) 

FOFX=XNCDF( X3,XMIJ,XSIG) 

IF( RN.LE.FGFX) go to 9125 
XA={ B+X3)/2. 

NBR=NBR+l 
CALL UNIF(RN,NU) 

F0FX = XNCDF( X4, X'lUtXSIG) 
IF(RN.LF.FOFX) GO TO 9063 
X5= { B+XA) /2. 

NBR=NBR+1 
CALL UNTP{RN,MU) 
FOFX=XNCDF( XEtXMU, XSIG) 
IFIRN.LE.FOFX) GO TO 1313 
EXMU=( B+X5) /2. 

EXSIG= ( X5-XA) /?. 

EXSIG=2. + '^XSIG 
£XSTG=EXSIG/6. 

GO TO 7000 

1313 EXMU={ X5+X4)/2. 

GXSIG={ X5-X^) 12 , 

EXSIG=2.+EXSIG 

6XSIG=EXSIG/6. 

GO TO 7C00 

9063 X5={ X3 + X2) /2. 

NBR=NBR+1 
CALL UNIF(RN,N'n 
FOFX=XNCDF( X5»XMU,XSIG) 
IF(RN.LE.FOFX) GO TO 1314 
X6=( X3<-X4) /2. 

NBR=NRR4-1 

CALL UNIF(RN,MU) 

FOFX=XNCDF( X6,XMU,XSIG) 

IF( RN.LE.FOFX) go to 1316 
EXMU=( X6+-X4) /2. 
EXSIG=(X6-X3)/2. 
exSTG=2.«EXSIG 
EXSIG==XSIG/6. 

GO TO 7000 
1316 EXMU=( X6+X3) /2. 

EXSIG=( X6-X3)/?. 

EXSIG=2 .*EXS1G 
EXSIG=FXSIG/6. 

GO TO 7000 

1314 X6=2.*X?-X5 
NBR=NBR+1 

CALL UNIF(PM,N'I) 
F0FX=XNCDF(X6,XWU,xSIG) 
IF(PN.LE.FOFX) GO TO 1315 
XB=X5 
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DFLX=X5-X2 

FXMU=XB+DELX/2. 

EXSIG=1 .3*DELX 
GO TO 7000 
1315 XB=X6 

nELX=X2-X6 
EXMU=XB-OFLX/'s-. 
FXSlG=iS>!<nELX 
GO TO 7000 
9125 X4= ( X2+X1 ) /2. 

N8R=NBR+1 
CALL UNIF(RN,MH) 
FOFX=XNCDF{ X^, XMIJ, XSIG) 
IF(PM.LE.FCFX) GO TO 909^ 
X5=( X2+X3) /2. 

NBR=NBR+1 

CALL UNIF(RNtNLI) 

FOFX=XNCDF( X5, X»^U,XSIG) 

IF( RN.LE.FOFX) GO TO 9047 
X6=( B4-X3) /2. 

NBR=NBR+1 
CALL UNIF(RN,NU) 
FOFX=XMCOF( X6,XMU,XSIG ) 
IF(RN.LF.FOFX) GO TO 9024 
X7=?.*B-X6 

nbr=nbr+i 

CALL UMIF(RO,Nlt) 

F0FX = XNICDF(X7, XMIJ , XSIG ) 

I F(PN.LE.FO>=^X) GO TO 9012 
XB=X7 
OELX-XB-B 
EXMU=XP+DELX/4. 
FXSIG=6.=)‘DELX 
GO TO 7000 
9012 XB=X3 

DELX=X6-XB 
EXf-MJ=XB + nELX/2. 
EXSIG=1.3*QELX 
GO TO 7000 
9024 EXMU=( X3+X5) /2. 

EXSIG=(X3-XF) f2, 

EXSIG = 2.’)'FXSIG 
6XSIG=FxSIG/6. 

GO TO 7000 
9047 X6={ X4+X2) /2. 

NeR=NBR+l 
CALL UNIF(RN,NIJ) 
FOFX=XNCDF( X6,XMU,XSIG) 
IF(PN.LE.FOPX) GO TO 9011 
FXMU=( X5+X2) /2. 

FXSIG=( X5-X2) /2. 
FXSIG=2.*FXSIG 
EXSIG=FXS IG/6. 

GO TO 7000 
9011 X7=2.*X4-X6 
N8R=NBR4-1 
CALL UMIF(RN,N0) 
FOFX=XNCDF( X7, XMU, XSIG) 
IF(RN.LE.FOFX) GO TO 9010 
XB = X6 

DELX=X?-XB 

FXMU=XB4-DELX/2. 

FXSIG= 1 .’-i'f^FLX 
GO TO 7000 
9010 XB=X7 

DELX=X4-X7 
EXMU=XB-OELX/A. 
EXSIG=6*0ELX 
GO TO 7000 
9094 X5=2*X1-X4 
NBR=N8R+1 
CALL UNIF(RN,NU) 
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FnFX = XNCDF( X5, XM(I, XSIG) 
IF(RN.LE.FOFX) GO TO 9003 
XB = X4 

DELX=X2-XA 
eXM0=XR+DFLX/2. 
EXSIG=r.3*DFLX 
GO TO 7000 

9003 XB=X‘i 
DELX=X1-X5 
EXMU=XB-DFLX/A. 
EXSIG=6.*DELX 
GO TO 7000 

9250 X3= { ;> + Xl )/2. 

NBR=NBR+l 

CALL UNIF(RN,NO) 

FOFX=XNCOF( X3,XM0,XSIG) 

IF( RN.LE.FOt^X) GO TO 93 75 
X4=( XI +X2 )/2. 

NBR=NBR+l 

CALL UNIF(RN,MU) 

FOFX = XMCOF( XA, XMU, XSIG ) 
IF(RM.LS.FOFX) GO TO 9004 
X5=(B+X2) /2. 

NBR=NBR+1 

CALL UNIFfRN.NU) 

FOFX = XNCDF( X5,XV(J,X9IG) 

IF( RN.LE.FOFx ) GO TO <^'005 

X6=2*B-X5 

NBR=NBR+1 

CALL UNIF(PN,NU) 

F0FX=XNCDF(X6,XMU,XSTG) 

IF(RN.LE.FO‘=X) GO TO 9006 

XB = X6 

DELX=X6-R 

EXMU=XB+nELX/4. 

EXSIG=6*DELX 

GO TO 7000 

9006 XR=X2 
DELX=X5-XB 
FXMU=XB+DELX/2 . 

EXSIG= l.3«DELX 
GO TO 7000 

9005 EXMU=( X2+X4) /2. 

EXSIG=( X2-X4)/2. 

FXSIG=2.^FXSTG 

EXSIG=EXSIG/6. 

GO TO 7000 

9004 X5=( Xl+X3)/2. 

NBR=NBR+1 

CALL UNIF(PM,NU) 

FOFX = X.NCDF( X5, XMU, X9IG ) 
TF(RN.LF.FOPX) GO TO 5555 
FXMU=( X1+X4) /2. 

EXSIG= (X4-X1 )/2. 

EXSIG=2.=<'FXSIG 

EXSIG=FXSIG/6. 

GO TO 7000 
5555 X6=2.*X3-X5 
NBR=NBR+1 
CALL ONIF(PN,NU) 
FOFX=XN'CnF( XA , XMO, XSIG ) 
IF(PN.LE.FOFX) GO TO 9007 
X8 = X5 

DELX = X l-xo 
EXMU=XB+nFLX/2. 

FXSIG=1 .3*DELX 
GO TO 7000 

9007 XB=X6 
DELX=X3-XB 
FXMU=XB-r>FLX/4. 
EXSTG=6*DELX 

GO TO 7000 
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9375 X4=?.*A-X3 
NBR=NBR+1 

CALL UNIF(PN,NU) 

FOFX = XNCDF( X4 , XMLI, XSI G ) 
IF(RN.LE.FOFX) GO TO 9376 
XB=X3 

OELX=Xl-XR 
FXMU=XR+OELX/2 . 

EXSIG=1 ,3*D£LX 
GO TO 7000 

9376 XB=X4 
DELX=A-XB 
eXMU=XR-DELX/4. 
EXSIG=6*DELX 

GO TO 7000 
9500 X2=(A+Xl)/2. 

NBR=NBR+l 
CALL UMIF(RNtNU) 
FOFX=XMCDF( X2 fXMU, XSIG) 
IF(RN.LE.FOFX) go TO 9501 
X3=( Xl+B) /?. 

NBP = NRR. + l 

CALL UNIF(RN»NU) 

FOFX=XNCCF( X3,XVU,XSIG ) 

IF(RN.LE.FOFX) GO TO 5556 

X4=2.*R-X3 

NBR=NBR+1 

CALL UMIF(RNtNU) 

FGFX = XNCDF( X4,XMIJ,XSIG ) 
IF(RN.LE.FOFX) GO TO 5554 
XR=X4 

DELX=X4-XR 
EXMU=XR + 0F.LX/4. 
EXSIG=6>^0ELX 
GO TO 7000 
5554 XR=Xl 

0ELX=X3-XB 
EXMU=XR+DFLX/2. 
EXSIG=1.3*DELX 
GO TO 7000 

5556 X4=( X1+X2) /2. 

NBR=NBR+1 

CALL UNIF(RN,NU) 
FOFX=XNCDF( X4,XMij,XSTG) 
IF(RN.LE.FOFX) GO TO 5557 
X5=( X1+X3) /2. 

N8R=NBR+1 

CALL UNTF(RN,NU) 

FOFX = X.NCDF( X5,XMt)»XSIG) 
IF(RN.LE.FO‘"X) GO TO 5553 
X6=2.*X3-X5 
NBR=NBR+l 
CALL UMIF(RM,MU) 
FOFX=XNCDF( X6, XMU, XSIG) 
IF(RN.LE.FOFX> GO TO 5559 
XR = X6 

DELX=XR-X3 
EXMU=XB+0ELX/4. 
EXSIG=6*DELX 
GO TO 7000 
5559 XB=X1 

DELX=X5-X1 
EXMU=XR+DELX/?. 
£XSIG=1.3+DELX 
GO TO 7000 
5553 EXMU=( X4+X1) /2. 

EXSTG= ( X]-X4) /2. 
EXSIG=FXS!G/3. 

GO TO 7000 

5557 X5=(A+X2)/2. 

NBR=NBR+l 

CALL UNIF(RN,NU) 
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FOFX=XMCDF( X5 , XMU , XSI G ) 
IP( RN.LE.FOFX) GO TO 5561 
X6=( X2+X4) /2. 

NBR=NBP+1 

CALL UNIF(RO,NH) 

FOFX = XNCDF( X6,XMI.|,XSIG ) 
IF(RM.LE.FOFXJ GO TO 121 
X7=( X^^ + X1 I /2. 

NBR=NBR+1 

CALL UNIF(RNJ,MU) 

FOFX=XNCD'^( X7, XMtJ, XSIG ) 

IFfRNI.LF.FOPX ) GO TO 122 

XR = 2 .-^Xl-X7 

nrr=nrr+i 

CALL UNIF(RN,^'U) 

FOFX = XNCDF(X8,X^’U,X5IG) 
IF(RN.LE.FOFXJ GO TO 123 
FXMU=XB+( XB-Xl > /A. 

EXSIG=6 .*( XS-Xl I 
GO TO 7000 

123 FXMl)= ( XA + X7) /2. 
EXSIG=1.3*(X7-X^^) 

GO TO 7000 

122 XR=( X6+X4) /2. 

NBP = N'BR + 1 
CALL UNIF(RN»Nll) 
FOFX=XNCOF( XS»XMM,XSTG ) 
IF(RN.LE.’"OFXJ GO TO 124 
X9=( X^+X7) /2. 

NRR=NPR+1 
CALL UNIP(RN,NU) 
FOFX=XNCOF( X9»XMU,XSIG) 
IF(RN.LE.FOFX) GO TO 125 
EXMLI=( X4 + X9) /2. 
EXSIG=1.3*(X9“X4) 

GO TO 7000 

125 FXMIJ=( X8+X4) /2. 

EXSIG=( X4-XR)/6. 

GO TO 7000 

124 X9=( X2+X6) /2. 
f^!8R=^3BR + l 

CALL UNIt^(RN,Nin 
FOFX=XNCOF(XQ, XMU,XSIG) 
IF(RN.LE.FOFXl GO TO 126 
EXMU={ X6+X8) /?. 

EXSIG=( X8-X6)/6. 

GO TO 7000 

126 EXMU=( X9^-X6) A2. 
EXSIG=1.?*( X6-X9) 

GO TO 7000 

121 X7=( X5+X2) /2. 

NRR=NBR+1 
CALL UNIF(RN,NU) 
F0FX=XNCDF(X7, XWU.XSIG) 
IF(RN.LF.FOFX) GO TO 127 
XB=( X6+X2) /2. 

NBR=MBR+1 

CALL UNIF(RN,N')) 

FOFX = XNC0F( X8 ,X'-OJ, XSIG) 
IF(RN.LE.FOFX) GO TO 128 
X9=( X6+X4) /2. 

NBR=NBR+-1 
CALL UN!F(RN,NII) 
FOFX=XNCDF( X9, XMIJ, XSIG) 
IF(RN.LF.FORX ) GO TO 129 
EXMU=( X6+X9) /2. 

FXSI G= 1 .3*( X9-X6) 

GO TO 7000 
129 EXMU=( XR+X6) /2. 

EXSIG=( X6-X8) /6. 

GO TO 7000 
128 X9=( X7+X2) /2. 
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NBR = NPR<-1 

CALL UN1F(RN,NII) 

FOFx = xNcnr( x9,XMti,xsir, ) 
IF(RN.LF.FOCX) on to 130 
FXMU=( X?fXB) /?. 

EXST G= ( XB-X2) /6. 

GO TO 7000 
130 EXMU=( X9+X2) /2. 

EXSIG=1.3*(X2-XO) 

GO TO 7000 

127 EXMI|=( X7 + X2) /2. 

EXSIG=1 .3«(X2-X7) 

GO TO 7000 
5561 X6=2.«A-X5 
NRR=NBR+1 
CALL U!MIF(PN,NU) 
FOFX=XNCDF( X6»XWI)t XSTG ) 
IF(RN.LE.FO>=X) GO TO 5560 
XB=X5 
DELX=X5-A 
EXMU=XB+0SLX/2. 
EXS!G=1.3«DELX 
GO TO 7000 
556C XB=X5 

D5LX=A-X6 
EXWU=XB-OELX/A. 
EXSIG=5*0ELX 
GO TO 7000 
9501 X3=(A+X2)/2. 

NRR=NBR+1 
CALL IJNIF(RM,NU) 
FOFX=XNCOF( X3»XMU,XSIG) 
IF(RN.LE.FOFX) GO TO 9503 
XA=( X1+X2) /2. 

NRR = N!BR + 1 
call 0NIF(RN,NU) 

FOFX = XNCDF{ XAtXMil, XSTG ) 
IF(RM.LE.FOFX) GO TO 9504 
X5=2 .^Xl-X4 

nbr=nrr+i 
CALL UNIF(RN,NU) 
FOFX=XNCDF( X5,XMU,XSIG) 
IF(PN.LE.FOFX) GO TO 9080 
XB = X5 

OELX=XB-Xl 
EXMU=XB+0ELX/4. 
EXSIG=6*0ELX 
GO TO 7000 
9030 XR=X2 

DELX=XA-X2 

FXMU=XR+0FLX/2. 

EXSIG=] .3’'<0ELX 
GO TO 7000 
9504 X5=( X3+X21/2. 

NBR=NRR+1 
CALL UNIF( RNJ ,Ni.n 
FOFX=XNCDF( X5»XMO,X5IG) 
IF(RN. LE.FOf^X) GO TO Q081 
X6={ X2+X4) /2. 

NBR=NBR+1 

CALL UNT‘^(RM,MO) 

F0FX = XNIC0F ( X6,XMU,XSIG) 
TF(RN.LE.FnFX) GO TO 9032 
X7 = 2 ,*X4-X6 
NBR=NBR+1 
CALL UNIF( RM,MO) 
FOFX=XNCDF (X7, X^O, XSIG) 
IF(RN.LE.FOFX) GO TO 9083 
XP=X7 

DELX=X7-X4 
EXMU = XR4-DELX/4. 
FX$IG=6*DFLX 
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Gn TO 7000 
9033 X3=X2 

DGLX=X6-X2 
FXMU=xR4-DELX/2 . 
EXSIG=1.R*DELX 
GO TO 7000 

9032 EXMU=( X5+X?) /2. 

FXSIG=( X2-XR) /6. 

GO TO 7000 
9031 X6=(X3+A)/2. 

NRR=N3R+1 
CALL UMIF(RM,MU) 
FOFX=XNCOF( X6,XVy»XSIG) 
IF(RN.LE.'=nFX) GO TO 903A- 
EXMU=( X3+X5) /2. 

EXSIG= ( X3-X3 ) /6. 

GO TO 7000 
903A X7=2.*A-X6 
MBR=NPR+l 
CALL UN)IF(RM,NU) 

F0FX = XMCDR(X7, XMl),X3TG ) 
IF(RN.LE.FOFX) GO TO 9035 
X3=X6 

DELX=X3-X6 

EXMU=XB+DELX/2. 

FXSTG= 1.3+DELX 
GO TO 7000 
0085 XB=X7 

0ELX=A-X7 
FXMU=XB-nFLX/^. 
EXSIG=6*D£LX 
GO TO 7000 
9503 X4=(X3+A>/2. 

NBR=NBR+1 
CALL UNIF{P^4t^)U) 

FOFX = XNCDF( X4,XMIJ, XSIG) 

IF (RN. LE.FOFX) GO TO 9507 
X5=( X2+X3) /2. 

NBR=NBR+l 

CALL UNTF(RO,Nt)) 

FOFX=XNCDF( X5.XMU, XSIG) 

IF(RN.LE.FOPX) GO TO 9509 

X6=2.*X2-X5 

r4BR = NBR + l 

CALL lJNIP(RN,Nin 

FOFX=XNCOF( X6,XMH,XSIG) 

IF(PN. LE.FOFX) GO TO 9510 

XB = X6 

D£LX=X6-X2 
FXMU=XB + nELX/A-. 
EXSIG=6*0ELX 
GO TO 7000 

9510 XB=X3 
DELX=X5-X3 
FXMU=XP+DELX/2. 

EXSIG=1 .3*0F.LX 
GO TO 7000 

9509 X6=( X3+X4) /2. 

NRR=MBR+1 
CALL UNIF(RN,NU) 

FOFX = XNCDF( X5tXM(),XSTG ) 
IF(RN. LE.FOFX) GO TO 9511 
FXMIJ=( X6 + X3) /? . 

FXSTG= ( X3-X6) /2. 

EXSIG=2.^'=XSIG 

EXSIG=EX5IG/6. 

GO TO TOOO 

9511 EXMU=( X44-X6) /?. 

EX3IG= ( X6-X4) /2. 

EXSIG=2.»EXSIG 

EXSIG='=XSIG/6. 

GO TO 7000 
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9507 X5=(X4+A)/2. 

NRR=NP04-1 

CALL UNIP(RM,NtJ) 

F0FX=XNCDF ( X5 , XMU, XSin ) 
IF(RN.LE.F0'=X ) on TO 9508 
EXMU=( X4+X5) /?. 

EXSIG= ( X4-X5) /2. 

EXSIG=2.*FXSIG 
EXSIG=EXSIG/6. 
on TO 7000 

9508 EXMU=( X5+A ) /2. 

EXS1G=( X4-X5) /2. 

EXSIG=2.4CXSIG 

EXSIG=EXSIG/6. 

GO TO 7000 

7000 X8=0. 

DELX=0. 

WRITE( 6, 100-i) FXMUtEXSIG 
IF(EXSIG.LT.O. ) EXSIG = -E.XSIG 
0PT=0. 

DIF=0. 

XINC=0. 

npT=EXMU-4.*EXSI G 
IF(OPT.LT.A) GO TO 1003 
XI NC=( nPT-A) /FXS IG 
INC=XINC/1 
DIF=XIMC-IMC 
IF(DIF.GE..5) INC=TNr+l 
IF(OIF .LT..5) INC=INC 
GO TO 1001 
1003 A=OPT 
INC = 0 

1001 M=I0 + IN'C + 1 
33 IS=70-N8R 

M=I0+INC+1 

CONDUCT 8RUCET0N TEST 

CLEAR ARRAYS 

DO 10 1=1,200 
X( I ) = 0. 

IXO( I )=0 
IXX( I )=0 
NS( I )=0 
NG( I ) = 0 
SUMAR=0. 

SUMBR=0. 

AR ( I )=0. 

BR( I )=0. 

10 CONTINUE 

LOAD X ARRAY 

DO 20 J=l,200 
X<J)=A + ( J-D+EXSIG 
20 CONTINUE 

CONDUCT experiment 
30 CALL UNIF(RN,nhi) 

FOFX=XNCOF( X( M) ,XMU,XSIG) 

IF(RN.GT.FOPX)GO TO 40 

IXX(M)=IXX(M)4-1 

M=M-1 

N=N+1 

50 IF(N.GT.IS) GO TO 60 
GO TO 30 

40 IX0(M)=IX0(M)<-1 

M = M+1 
N=N+1 
GO TO 50 



PERFORM 8RUCET0N ANALYSIS 
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COUNT RESPONSES AMD NON-RESPONSES 
6)0 ISIIMX=0 
ISUNn=o 
on lA J=l,200 
ISUMX=ISUMX<-IXX( J) 

ISUMO=I SUMO+IXO( J) 

NS( J ) =0 
AR ( J ) = 0 . 

RR( J ) = 0. 

NG( J )=0 

lA CONTINUE 

DETERMINE LESS FREQUENT EVENT AND LOAD NS 
IF( isumx.le.isumq) go to 15 
NT=I SUMO 
IFLAG=0 
DO 21 J=lr200 
NS( J)=IXO( J) 

21 continue 

GO TO 16 
15 NT=ISUMX 
IFLAG=1 
DO 2 2 J= 1,200 
NS( J ) = IXX( J) 

22 CONTINUE 

DETERMINE FIRST AND SECOND MOMENTS 
If JC0UNT=1 

17 IF{ NS( vICOUNT) . GT.O) GO TO IS 
JCOUNT=JCOUNT+l 
IFUC0UNT.GF.2O0) GO TO 10«^ 

GO TO 17 

IB MC0UNT = 200-.IC0iiNT 
DO 19 J=1,MC0UMT 
NG( J)=NS( JC0IJNT4J-1) 

AR (J) = U-1 I’t'NGU ) 

SUMAR=SUMAR+AR{ J ) 

BR(J ) = ( ( J-1 )’?^*2)*NG( J) 

SUMPR=SUMPR+BR ( J ) 

19 CONTINUE 

YPRIME = XUCOUNT) 

CALCULATE ESTIMATES 0'= MEAN AND STANDARD DEVIATION 

IF ( IFL AG.FO.O) XMUEST = YPPI MF+EXSI G*{ ( SUMAR/NT ) 4-( 1./2. ) ) 
IF ( .NOT. IFLAG.EO.O) XMUEST=YPRI MF+exSIG*( ( SUMAR/NT )-( 1 . 
SIGFAC=( (NT*SUMRR)-(SUMAR + -i'2) ) / (NT*^*2) 

IFISIGFAC.GT.. 3) GO to 1000 
FMU(LC0UNT)=0. 

EOEV(LCOUMT)=0. 

N0G0=N0G04-1 
GO TO 104 

1000 DEVEST=l.f2’'"EXSIG*( SIGFAC+.029) 

LOAD EMU AND EDEV 

FMU( LCOUNT)=XMUFST-XMU 
EDEV(LCOUNT) =DFVEST-XSIG 
AODMU=ADDMU+EMU( LCOUMT) 

ADDSIG=ADDSIG+EDEV{ LCOUNT) 

ADDMU0=ADDVU04-FMU( LCOUNT ) v* 2 
ADDS DO = ADDSDO+FDEV( LCOUNT) 1^*2 
IF(FMU(LC0UMT) .LT. 0. ) GO TO 91 
IF(EMU(LCnUNT) .EO.O. ) GO TO 92 
IMUHI=IMUHI+1 

IF(FmU(LCOUMT) .GT.HIMU) HI Mij=F MU ( LCOUMT ) 

IF( .NOT.EMU(LCOUNT) .GT.HIMU) HI MU=HImij 
GO TO 93 

92 N0MU=N0MU4-1 
GO TO 93 

91 1MUL0=IMUL04-1 

I F(SMU( LCOUNT) .LT.SMLO) SMLO=E MU (L COUNT ) 
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IF( .NOT.EMU( LCOUNT) .LT. SMl.ni SMl.O=SMLO 
^>3 IF(FDEVaCOUNT) .IT.O. ) GO TO 9 ^ 

IF{ FDEV( LCOUNT) .eO.O. ) GO TO 05 
IDEVHI=lnFVHl+l 

IF ( EOF V( LCOUNT) .GT.OEVH] ) DFVH I =FDFV ( LCOUNT ) 

IF( . NOT. EOEV( LCOUNT) .GT.DFVHI ) nEVHI=DEVHI 
GO TO 104 

95 NDDEV=N00EV4-1 
GO TO 10^ 

94 IDEVLO= IDFVLO + 1 

IFU:DEV(LC0')NT) .LT.OEVLO) DEVLO=EOEV ( LCOUNT ) 

IF( . NOT. EOEVI LCOUNT ). LT.OEVLO) OEVLO=OFVLO 
104 JCOUNT=0 
SIGFAC=0. 

XMUFST=0. 

DEVF5T=0. 

SUMAR=0. 

SU^'BR = 0. 

LC0UNT=LC0UNT^1 

HAVE 1000 EXPERIMFMTS SEEN CONDUCTED ? 
IFILCOUNT.LT.lOO] ) GO TO 103 
IF EXPERIMENTS COMPLETED CALCULATE AND WRITE P^^SULTS 
EXNOGO=NOGO 

SAMAVM=AODMU/( 1 000. -F XNOGO ) 

SAMAVD=A0DS IG/( 1000. -F XNOGO) 

SAMS OM = AODMUO/ ( QT9.-F XNOGO) 

SAMS00=ADD$n0/(99<3 .-F XNOGO) 

IFUANY.EO.l) GO TO 35 
IANY=I ANY+l 
GO TO 69 
35 STOP 
END 



SUBROUTINE UNIF{RN,NU) 



SUBROUTINF RETURNS RANDOM NUMBER UNIFORM ON (0,1). 

real MOD 
MOD= 2<'*3l 
NR=129*NU+1 
RN=NR/MOD 

IF(RN.LT.O.O) rn=-rn 

NU=NR 

RETURN 

END 



FUNCTION XNCDF( V,XMU,SX) 



FUNCTION SUSPRTGRAM CALCULATES CUMULATIVE NORMAL. 

X IS AN R.V. WITH MEAN, XMU, AND STANDARD DEV I AT I ON , S X . 
ARG=( V-XMU)/SX 
XNCDF=SNCDF( ARG ) 

RETURN 

END 



function SNCDF(X) 
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FUNCTION SUBPROGRAM CALCULATES STANDARD CUMULATIVE NORMAL. 
DATA TFST/0.0/ 

IF( TEST. NE .0.0) GO TO 100 
SR2= SQRT(2.0) 

T5ST=1 . 

100 SNCDF=(1.0+FRF(X/SR2))/2.0 

R'^.TURN 
END 
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